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We theoretically investigate the phonon scattering by vacancies, including the impacts of missing mass and 
linkages (xy 1 ) and the variation of the force constant of bonds associated with vacancies (x^ 1 ) by the 
bond-order-length-strength correlation mechanism. We find that in bulk crystals, the phonon scattering rate 
due to change of force constant x^ 1 is about three orders of magnitude lower than that due to missing mass and 
linkages Xy 1 . In contrast to the negligible x^ 1 in bulk materials, x^ 1 in two-dimensional materials can be 3-10 
folds larger than Xy 1 . Incorporating this phonon scattering mechanism to the Boltzmann transport equation 
derives that the thermal conductivity of vacancy defective graphene is severely reduced even for very low 
vacancy density. High-frequency phonon contribution to thermal conductivity reduces substantially. Our 
findings are helpful not only to understand the severe suppression of thermal conductivity by vacancies, but 
also to manipulate thermal conductivity in two-dimensional materials by phononic engineering. 



Thermal property of two-dimensional materials has attracted considerable attention due to the unique 
physical attributes 1 " 4 and the potential applications 5 " 8 in many areas. On the one side, because of its superior 
high thermal conductivity, graphene-based devices can be applied in nanoscale phononics and for thermal 
management 9 " 16 . On the other side, the ultra low thermal conductivity of monolayer molybdenum disulfide 
(MoS 2 ) sheet and nanoribbon 17 " 19 opens up the possibility to realize two-dimensional thermoelectric devices 
based on Transition metal dichalcogenides (TMD) materials for heat energy harvesting and refrigeration. Most 
materials have natural point defects, such as atomic vacancies introduced in the process of fabrication. Since the 
phonon frequency depends on length and energy of the local bonds 20 , these vacancy defects offer a possibility for 
tailoring the thermal conductivity of materials 21 " 25 . 

Phonon scattering by vacancies in crystals was originally treated from perturbation theory by Ratsifaritana and 
Klemens. Although the role of the strain field in scattering of phonons is important since strain would change the 
local value of phonon frequency for a fixed wave vector, thus contributing to the perturbation, Ratsifaritana and 
Klemens demonstrated that the distortion effects could be neglected for vacancies, and the perturbation of 
vacancy defect was one that corresponds to the removal of the mass of one atom and the force constants of 
two atoms 26 . In addition to the missing mass and missing linkages (Ty 1 ), the change of force constant of bonds 
between the under-coordinated atoms near the vacancies also results in phonon scattering (t^ 1 )' t> ut & is 
neglected in the theoretical model of vacancy. Usually, the interatomic force constants can be obtained by the 
quantum chemistry calculations. However, it is not reliable to investigate the vacancy effect on phonon scattering 
in sizeable materials due to computational limitation. 

In this communication, we report our findings on the phonons scattering (both Xy 1 and x^ 1 ) by vacancies, 
including the impacts of missing mass, missing linkages, and change of force constant of the under- coordinated 
atoms near the vacancies based on the frame work of bond- order-length- strength (BOLS) notation 27,28 . We find 
that the phonon scattering rate due to change of force constant x^ 1 is about three orders of magnitude less than 
that due to missing mass and linkages Xy 1 in bulk crystals. In contrast to the negligible x^ 1 in the bulk, x^ 1 in two- 
dimensional materials can be 3-10 folds larger than Xy 1 . The critically different phonon scattering rate deter- 
mines that in two-dimensional materials with vacancies, the contribution to thermal conductivity from high 
frequency phonons is greatly suppressed. Consistency in the present and molecular dynamics simulations 29 on 
the thermal conductivity of defected single-layer graphene confirmed our expectations. 
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Figure 1 | Schematic description of graphene flakes with (a) single vacancy. The CN of carbon atoms in graphene is 3. There are 3 under-coordinated 
atoms displayed as red colour whose CN is 2 near the single vacancy, and bonds displayed as blue colour between the under-coordinated atoms near the 
vacancy become stronger, (b) double vacancy. 



Results 

Bond-order theory for phonon-vacancy scattering. Perturbation 
theory in terms of the missing mass and the missing linkages deals 
with the scattering of phonons by vacancy defects in crystals 26 ' 30 , 



-i ( AM 2 n(o 2 g(co) 

T v = X( ) 

v M 2 G 



(1) 



where x is the density of vacancies, G is the number of atoms in the 
crystal, g(co) is the phonon density of states (DOS). For vacancy 

AM M a 

defect, the effective value of is — 2, where M is the 

M M 

average mass per atom, M a is the mass of the missing atom, and 
the term —2 accounts for the potential energy of the missing 
linkages, or twice the potential energy per atom 26,31 . 

One important fact that has been overlooked in previous modeling 
is that bond between under- coordinated atoms becomes shorter and 
stronger 28 . Bond shortening and strengthening not only raises the 
local density of charge, mass and energy but also deepens the local 
potential, providing perturbation to the local potential. This under- 
coordination effect raises the local energy density, or the elastic 
modulus, and subsequently the trap for phonon and electron trans- 
portation 32,33 . The coordination number (CN) dependence of bond- 
contraction coefficient (C z ) and bond energy (E z ) follow the 
relation 28 , 

{C z = d z /do =2/{l + exp[(12 — z)/(Sz)]} (bond — contraction coefficient) 
E z = C~ m E b (Single - bond - energy) ( 2 ) 

where z is the effective coordination number, d z is the bond length. In 
BOLS theory, the effective CN for bulk materials is always 12, regard- 
less of the nature of the bond or the crystal configuration 34 , E h and d 0 
are the corresponding bulk values of single bond energy and bond 
length, respectively, m is an indicator for bond nature of a specific 
material, which is not freely adjustable. For alloys and compounds 
(such as MoS 2 and BN) m is around 4, for C and Si the value of m has 
been optimized to be 2.56 and 4.88, respectively 35 . 

Eq. (2) yields the change of force constant of bonds between 
under-coordinated atoms, from the perspective of dimensionality 
analysis 28 , 



k z = d 2 u(r) /dr 2 \ r=dz ccE z /d 2 z 



(3) 



Accordingly, the bond force constant nearby the vacancies increases 
compared with that far away from vacancies. As shown in Fig. 1 (the 
schematic description of graphene flakes with single and double 
vacancy), the bonds displayed as blue colour between the under- 
coordinated atoms near the vacancy become stronger. According 



to the perturbation theory of Klemens, the rate of phonon scattering 
by atoms of different force constant is given as 36 ' 37 , 



,dk^ w 2 g((jo) 



(4) 



where x is the density of imperfections, k is the force constant, and 5k 
is the change of force constant. 

According to Eq. (2) and Eq. (3), we can get 

k z -i = £z 
k z 

(5) 
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where k z - 1 is the force constant of under-coordinated atoms near the 
vacancies whose CN is z-l, and k z is the force constant of atom whose 
CN is z. Therefore, according to Eq. (4) and Eq. (5), the scattering rate 
of phonons by the under- coordinated atoms near the vacancies is 
derived as, 
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where X A — zx IS the density of the under-coordinated atoms, and x is 
the density of single vacancy, m is a key parameter that represents the 
nature of the bond. 

The unusual vacancy effects on phonon scattering in 2D materials. 

The ratio of t^ 1 to Ty 1 in the case of different m and different 
coordination number is shown in Fig. 2. For bulk materials whose 
effective atomic CN is always 12 in the BOLS correlation mechanism, 
t^ 1 is three orders of magnitude less than t^ 1 , so the scattering of 
phonons by the under-coordinated atoms near the vacancies is 
negligible in bulk materials. The ratio t^ 1 j%y X increases quickly 
with decreasing the coordination number. For two-dimensional 
materials with z = 3, such as silicene, hexagonal Boron Nitride, 
graphene and MoS 2 , t^ 1 is 3- to 10-fold larger than Ty \ therefore 
the scattering of phonons by the under- coordinated atoms near the 
vacancies must be taken into account, and the total scattering rate of 
phonons by the imperfections should be Ty 1 + tT~ 1 . 

Moreover, the ratio t A ~ 1 / T v~ 1 increases with the parameter m 
increases. For carbon-based 2D materials, that is, graphene, m is 
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Figure 2 | The ratio of 1 to Xy 1 in the case of different m and different 

CN. Coordination number z=\2 represents the phonon scattering in bulk 
materials. 
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Figure 3 | The phonon scattering rate as a function of phonon angular 
frequency. The phonon scattering mechanisms including Umklapp, 
boundary, vacancy and under-coordinated atoms scattering are considered. 



found to be 2.56. From Fig. 2, the ratio 1 /Ty 1 is 2.8. For 2D alloys 
and compounds, such as BN and transition metal dichalcogenides 
(TMDs) materials, m is around 4, results to a high t^ 1 /iy 1 ratio of 
6.4. Very recently, monolayer MoS 2 , a member of TMD family, have 
gained considerable interest. Thus for 2D TMDs materials, the pho- 
non scattering by under- coordinated atoms is dominated over that 
by vacancy itself. 

Vacancy effects on phonon scattering in graphene flakes. It is an 

important phononic engineering technique to modulate the thermal 
transport by phonons with different range of frequency 38 . Next using 
graphene flake as an example, we discuss the impacts of vacancies on 
phonon scattering rate and thermal conductivity contributed from 
phonons with different frequency. The Matthiessen's rule which 
assumes that different scattering mechanisms are independent is 
adopted, so the total phonon scattering rate 1 in branch X is given as: 

(7) 



l u\ + T B,l + T V,\ + T A,l 



The Umklapp phonon-phonon scattering rate (i u \) and the phonon 
boundary scattering rate ) are given as 39 , 

ylk B Tw 2 



Mv\cd d ,x 



W 1+P 



(8) 



(9) 



where M is the mass of a graphene unit cell, yx is the Griineissen 
parameter, W is the width of graphene ribbon, is the component 
of the phonon velocity in branch X perpendicular to the longitudinal 
direction of graphene ribbon, and P is the specularity parameter, which 
is defined as the probability of phonon's specular reflection at the lateral 
boundaries. In order to fit the experimental thermal conductivity value 
of pristine graphene, we set the specular parameter as 0.8, and the width 
of graphene flake as 5 |im, which are consistent with Ref. 40 and Ref. 41. 
As the major concern of the present work are the mechanism of 
phonon scattering by vacancies and the reduction of thermal 
conductivity by vacancies, the values of specular parameter and flake 
width are fixed. The Griineissen parameters for acoustic branches 
originate from the results of first-principles calculations in Ref. 42. 
Based on the linear dispersion for in-plane branches and quadratic 
dispersion for out-of-plane branch (see Eq. 16 in the section of 
Methods), the function for phonon DOS of graphene is, 



g(co)=< 



( Sea 

2nvl 
S 

I 4not 



X = LA, TA 
1 = ZA 



(10) 



Therefore according to Eq. (1), the rate of phonon scattering by single 
vacancies is given as, 



l v\~- 



2.25xQ- 



1.125x0- 



X = LA, TA 
k = ZA 



(11) 



where x is the density of vacancies, and Q is the primitive cell area of 
graphene. 

As shown in Fig. 1(a), the CN of carbon atoms in graphene is 3. 
There are 3 under-coordinated atoms whose CN is 2 near the single 
vacancy. Based on BOLS theory, the bond between the under-coordi- 
nated atoms becomes shorter and stronger. Girit et al 43 discovered that 
breaking a C-C bond of the 2-coordinated carbon atoms near the 
monolayer-GNR vacancy required 7.50 eV per bond, that was 32% 
higher than the energy (5.67 eV/bond) required for breaking one bond 
between 3-coordinated carbon atoms in the interior of a suspended 
graphene sheet. The mechanical strength of graphene increases with 
the density of defects 44 (the reconstructed 5- and 7-atom rings forming 
the grain boundaries), because of the particular strength of the ring 
bonds and their elongation dynamics. These findings provide evidence 
for the BOLS prediction of the shorter and stronger bonds at vacancies. 
According to Eq. (5), the ratio of force constant of the 2-coordinated 
carbon atoms near the vacancy to that of the 3-coordinated carbon 
atoms is 2.03, and according to Eq. (6), the rate of phonon scattering by 
the 2-coordinated carbon atoms near the single vacancies is given as, 



6.36x0- 



3.18x0- 



X = LA,TA 
X = ZA 



(12) 



From Eq. (12) and Eq. (1 1), we know that the phonon scattering by the 
under-coordinated atoms near the vacancy dominates over that by the 
vacancy in graphene. 

In the pristine graphene sheet, the Umklapp phonon-phonon 
scattering is dominant at room temperature, but in the defective 
graphene, the situation is quite different. Fig. 3 presents the fre- 
quency-dependent scattering rate of LA phonons for all kinds of 
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Figure 4 | The ratio of the thermal conductivity of defective graphene to 
the pristine one at room temperature as a function of defect density. The 

squares are the molecular dynamics simulation results of single vacancy 
defective graphene from Ref. 29. 

mechanisms in the single vacancy defective graphene. The vacancy 
density is 1%. For the ultra-low frequency phonons, the boundary 
scattering is dominant. However, for w > 2.5 Trad/s, the scattering 
by the under- coordinated carbon atoms dominates over other 
mechanisms. The situation is similar for other branches. However, 
this most important scattering mechanism is neglected in previous 
theoretical model on phonon scattering by vacancies. 

Vacancy effects on thermal conductivity in graphene flakes. Many 
researches presented that the thermal conductivity of graphene could 
be remarkably reduced by point defects 22 ' 29,45 , such as isotope and 
vacancy. The vacancy effect on the thermal conductivity is more 
remarkable than that of isotope. By molecular dynamics simula- 
tion, Haskins et at. found that even for a very low vacancy 
concentration 0.1%, a 81% reduction of thermal conductivity of 
graphene was achieved 29 . Although Klemens's model agrees well 
with thermal conductivity reductions by vacancies in bulk 
materials 31,46 , the severe suppression of thermal conductivity of 
graphene by vacancies is much over the prediction of Klemens's 
theory 22 . Incorporating phonon scattering by the under-coordi- 
nated atoms near the vacancies to the linearized phonon 
Boltzmann transport equation within relaxation time approxima- 
tion (see in the section of Methods), we calculate the relative 
thermal conductivity (k/k 0 ) of defective graphene (single vacancies 
and double vacancies) as a function of defect density at room 
temperature, which is shown in Fig. 4. k 0 (the thermal conducti- 
vity of pristine graphene) calculated by this work, 3750 W/mK, 
is coherent with various experimental measurements 47,48 and 
theoretical calculations 49,50 . The dash line is the calculation results 
by Klemens's model from Ref. 26, in which the phonon scattering by 
the under-coordinated carbon atoms near the vacancies is not 
considered. The solid lines are the results of the present model. 
The squares are the molecular dynamics simulation results of 
single vacancies in graphene from Ref. 29. It is obvious that the 
results predicted by the present model agrees better with molecular 
dynamics simulations than the results of Klemens's theory, and the 
phonon scatterings by the 2-coordinated carbon atoms near the 
vacancies further suppress the thermal conductivity of graphene. 

Fig. 5 demonstrates the normalized accumulative distribution of 
thermal conductivity defined as Eq. (19) (see in the section of 
Methods) for different single vacancy density. The curves move towards 
the upper-left corner as the vacancy density increases. The ultra-low 
thermal conductivity of vacancy defective graphene stems from the 
significant suppression of high-frequency phonons. With the increasing 
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Figure 5 | The normalized accumulative distribution of thermal 
conductivity of graphene flakes as a function of phonon angular 
frequency. Here the width of graphene flake is 5 (am, concentrations of 
single-vacancy defect are 0%, 0.01%, 0.1%, and 1%, respectively. The solid 
lines are results from the present model. The dotted line is thermal 
conductivity of graphene flake with 1% vacancy calculated without 
considering the phonon scattering by under-coordinated atoms. 

of vacancy density, the contribution of more and more high-frequency 
phonons to heat conductance is severely suppressed by the scattering of 
the vacancies and the under-coordinated atoms near the vacancies, as 
shown in Fig. 5. From Fig. 3, we know that the scattering by these 
imperfections is dominant, except for the ultra-low frequency phonons. 
Therefore, the thermal conductivity of vacancy defective graphene is 
absolutely dominated by low-frequency acoustic phonons. For example, 
in the case of 1% concentration of single vacancy, the phonons of co < 
3 Trad/s have about 50% contribution to the thermal conductivity, 
which is shown in Fig. 5. However, in the Klemens's model (the dotted 
line), the same phonons only have about 30% contribution which is 
lower than our model, because it neglects additional scattering by the 
under-coordinated atoms near the vacancies. 

Discussion 

The theoretical model in present work can be used to explain the 
different content dependence in thermal conductivity of graphene 
flakes with single- and double-vacancy defects. For a single vacancy, 
there are 3 under- coordinated atoms around each defective atom, 
while there are only 2 under- coordinated atoms around each defect- 
ive atom for a double vacancy, as shown in Fig. 1(b). According to Eq. 
(6), the rate of phonon scattering by the 2-coordinated carbon atoms 
near the double vacancies at the same defect density x is given as, 

{4.24x0^- 1 = LAJA 
2 (13) 
2.12x0— k = ZA 
a 

where the defect density x is defined as the number of defected atoms 
divided by the total atom number in pristine graphene. The rate of 
phonon scattering by the 2-coordinated carbon atoms near the single 
vacancies is as 1.5 times as that near the double vacancies at the same 
defect density. Therefore, the single vacancy defective graphene pos- 
sesses the lower thermal conductivity than the double vacancy defective 
one. For example, for graphene flakes with 0.1% defected atoms, for 
single vacancy defective graphene, k/k 0 is 0.19, while for double 
vacancy defective graphene, k/k 0 is 0.23. Thus under the same defected 
atom density, single-vacancy can result in larger reduction in thermal 
conductivity with respect to double-vacancy defects. This explains the 
phenomenon observed by molecular dynamics simulations 29,51 . 
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In summary, we find a new phonon scattering mechanism which 
is the scattering by the under-coordinated atoms near the vacancies. 
This mechanism origins from the increase of force constant of bonds 
associated with vacancies. The scattering of phonons by the under- 
coordinated atoms near the vacancies is negligible in bulk materials. 
However, because of the low coordination number, this mechanism 
has dominant effect on the phonon transport in two-dimensional 
vacancy defective materials. This finding is helpful not only to under- 
stand the severe suppression of thermal conductivity by vacancies, 
but also to apply phononic engineering to manipulate thermal con- 
ductivity in two-dimensional materials. 

Methods 

According to the linearized phonon Boltzmann transport equation within relaxation 
time approximation, the thermal conductivity in branch X of single layer graphene 
(SLG) in the y direction (the longitudinal direction of graphene ribbon) is derived as: 



Kx= -^f\ c t hV ly XiA ^ 



(14) 



where S is the area of the sample, vx, y is the y component of the group-velocity vector 
in branch X (1 = LA, TA and ZA, only acoustic branches are considered 49 ' 52 ), Xx is the 
averaged phonon relaxation time between successive scattering events, q is the wave 
vector, and c ph is the volumetric specific heat of each mode, which is given as: 



k B (Hco/k B T) 2 e nc °/ kBT 



(15) 



where k B is the Boltzmann constant, S = 0.335 nm is the thickness of graphene, H is 
the reduced Planck constant, and T is the absolute temperature. 

In SLG, the LA and TA acoustic branches are linear, whereas the ZA branch shows 
a quadratic dependence of the frequency on the wave vector 49,53 , so 



aq A 

Using this dispersion and the relationship vx - 
wave vector (q) to frequency (co), we can get 

k B [\hco/k B T) 2 e Hco / ksT 



X = LA, TA 
X = ZA 



(16) 



dq 



, and transforming the integral of 
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where C0d,a is the Debye frequency, which is given as 



(i)X)d(j) X = ZA 



2vx 
47ia 



X = LA,TA 



l = ZA 



(18) 



where Q is the primitive cell area. 

The total phonon scattering rate xj x is the sum of x^\, x^j, Xy\ and t^~], which 
are given as Eq. (8-12) respectively. For a flat graphene sheet lying in the x-y plane, the 
reflection symmetry requires that the Hamiltonian be invariant under z — » — z 54 . Seol 
et al. obtained a selection rule for three-phonon scattering, which requires that an 
even number of ZA phonons is involved in each process. This selection rule is adopted 
in our calculation. 

The normalized accumulative distribution of thermal conductivity is defined as, 



F(cd)-- 



K LA {CJ) + K T a(co) + K Z a(m) 
KlA + Kta + KzA 



where Kx{co) is given as, 



Kx(oS) - 
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